home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Language/OS - Multiplatform Resource Library
/
LANGUAGE OS.iso
/
cpp_libs
/
rwvector.lha
/
RWVector2.1
/
src
/
xvecmath.cc
< prev
next >
Wrap
C/C++ Source or Header
|
1989-08-18
|
7KB
|
304 lines
/*
* Math definitions
*
* Copyright (C) 1988, 1989.
*
* Dr. Thomas Keffer
* Rogue Wave Associates
* P.O. Box 85341
* Seattle WA 98145-1341
*
* Permission to use, copy, modify, and distribute this
* software and its documentation for any purpose and
* without fee is hereby granted, provided that the
* above copyright notice appear in all copies and that
* both that copyright notice and this permission notice
* appear in supporting documentation.
*
* This software is provided "as is" without any
* expressed or implied warranty.
*
*
* @(#)xvecmath.cc 2.1 8/18/89
*/
#include "rw/<T>Vec.h"
#define TYPE <T>_TYPE
#include "vecdefs.h"
// The type returned by abs depends on the type of vector
#if TYPE==DComplex_TYPE || TYPE==FComplex_TYPE
#define ABS_TYPE <P>
#define ABS_VEC <P>Vec
#else
#define ABS_TYPE <T>
#define ABS_VEC <T>Vec
#endif
#ifdef APPLY_MATH
#if TYPE==DComplex_TYPE || TYPE==FComplex_TYPE
<T>Vec
<T>Vec::apply(CmathFunTy f)
// Apply function returning <T> to each element of vector slice
{
register i = length();
<T>Vec temp(i);
register <T>* sp = data();
register <T>* dp = temp.data();
register j = stride();
while (i--) { *dp++ = f(*sp); sp += j; }
return temp;
}
<P>Vec
<T>Vec::apply2(CmathFunTy2 f)
// Apply function returning <P> to each element of vector slice
{
register i = length();
<P>Vec temp(i);
register <T>* sp = data();
register <P>* dp = temp.data();
register j = stride();
while (i--) { *dp++ = f(*sp); sp += j; }
return temp;
}
// Version 2.0 has trouble finding the proper overloaded real & imag.
// Put in an explicit version:
<P>Vec
real(const <T>Vec& v)
{
register i = v.length();
<P>Vec temp(i);
register <T>* sp = v.data();
register <P>* dp = temp.data();
register j = v.stride();
while(i--) {*dp++ = real(*sp); sp += j;}
return temp;
}
<P>Vec
imag(const <T>Vec& v)
{
register i = v.length();
<P>Vec temp(i);
register <T>* sp = v.data();
register <P>* dp = temp.data();
register j = v.stride();
while(i--) {*dp++ = imag(*sp); sp += j;}
return temp;
}
#else // Not Complex
<T>Vec
<T>Vec::apply(mathFunTy f)
// Apply function returning <T> to each element of vector slice
{
register i = length();
<T>Vec temp(i);
register <T>* sp = data();
register <T>* dp = temp.data();
register j = stride();
while (i--) { *dp++ = f(*sp); sp += j; }
return temp;
}
#endif
#endif
ABS_VEC
abs(const <T>Vec& s)
// Absolute value of <T>Vec slice.
{
register i = s.length();
ABS_VEC temp(i);
register <T>* sp = s.data();
register ABS_TYPE* dp = temp.data();
register j = s.stride();
while (i--) { *dp++ = ABSOLUTE(*sp); sp += j; }
return temp;
}
#if HAS_ATAN2
<T>Vec
atan2(const <T>Vec& u,const <T>Vec& v)
// Arctangent of u/v.
{
register i = v.length();
u.lengthCheck(i);
<T>Vec temp(i);
register <T>* up = u.data();
register <T>* vp = v.data();
register <T>* dp = temp.data();
register uj = u.stride();
register vj = v.stride();
while (i--) { *dp++ = atan2(*up,*vp); up += uj; vp += vj; }
return temp;
}
#endif
<T>Vec
cumsum(const <T>Vec& s)
// Cumulative sum of <T>Vec slice.
// Note: V == delta(cumsum(V)) == cumsum(delta(V))
{
unsigned i = s.length();
<T>Vec temp(i);
register <T>* sp = s.data();
register <T>* dp = temp.data();
REGISTER <T> c = 0;
register j = s.stride();
while (i--) { c += *sp; *dp++ = c; sp += j; }
return temp;
}
<T>Vec
delta(const <T>Vec& s)
// Element differences of <T>Vec slice.
// Note: V == delta(cumsum(V)) == cumsum(delta(V))
{
unsigned i = s.length();
<T>Vec temp(i);
register <T>* sp = s.data();
register <T>* dp = temp.data();
REGISTER <T> c = 0;
register j = s.stride();
while (i--) { *dp++ = *sp - c; c = *sp; sp += j; }
return temp;
}
<T>
dot(const <T>Vec& u,const <T>Vec& v)
// Vector dot product.
{
register i = v.length();
u.lengthCheck(i);
REGISTER <T> t =0;
register <T>* up = u.data();
register <T>* vp = v.data();
register uj = u.stride();
register vj = v.stride();
while (i--) { t += *up * *vp; up += uj; vp += vj; }
return t;
}
#if HAS_MINMAX
int
max(const <T>Vec& s)
// Index of maximum element.
{
if (s.length() == 0) s.emptyErr("max");
register <T>* sp = s.data();
REGISTER <T> t = *sp;
register j = s.stride();
register x = 0;
for (register i =0; i<s.length(); i++) {
if (*sp > t) { t = *sp; x = i; }
sp += j;
}
return x;
}
int
min(const <T>Vec& s)
// Index of minimum element.
{
if (s.length() == 0) s.emptyErr("min");
register <T>* sp = s.data();
REGISTER <T> t = *sp;
register j = s.stride();
register x = 0;
for (register i =0; i<s.length(); i++) {
if (*sp < t) { t = *sp; x = i; }
sp += j;
}
return x;
}
#endif
<T>
prod(const <T>Vec& s)
// Product of <T>Vec elements.
{
register i = s.length();
register <T>* sp = s.data();
#if TYPE==DComplex_TYPE || TYPE==FComplex_TYPE
REGISTER <T> t = <T>(1,1);
#else
REGISTER <T> t = 1;
#endif
register j = s.stride();
while (i--) { t *= *sp; sp += j; }
return t;
}
#if HAS_POW
<T>Vec
pow(const <T>Vec& u,const <T>Vec& v)
// u to the v power.
{
register i = v.length();
u.lengthCheck(i);
<T>Vec temp(i);
register <T>* up = u.data();
register <T>* vp = v.data();
register <T>* dp = temp.data();
register uj = u.stride();
register vj = v.stride();
while (i--) { *dp++ = pow(*up,*vp); up += uj; vp += vj; }
return temp;
}
#endif
<T>Vec
reverse(const <T>Vec& s)
// Reverse vector elements.
{
register i = s.length();
<T>Vec temp(i);
register <T>* sp = s.data();
register <T>* dp = &temp(i-1);
register j = s.stride();
while (i--) { *(dp--) = *sp; sp += j; }
return temp;
}
<T>
sum(const <T>Vec& s)
// Sum of a <T>Vec
{
register int i = s.length();
register <T>* sp = s.data();
REGISTER <T> t = 0;
register j = s.stride();
while (i--) { t += *sp; sp += j; }
return t;
}
#if HAS_VARIANCE
<P>
variance(const <T>Vec& s)
// Variance of a <T>Vec
{
register int i = s.length();
if(i<2){
RWnote("::variance()", "Less than two points.");
RWerror(DEFAULT);
}
register <T>* sp = s.data();
register <P> sum2 = 0;
REGISTER <T> avg = mean(s);
register j = s.stride();
while (i--) {
REGISTER <T> temp = *sp - avg;
#if TYPE==DComplex_TYPE || TYPE==FComplex_TYPE
sum2 += norm(temp);
#else
sum2 += temp * temp;
#endif
sp += j;
}
// Use the biased estimate (easier to work with):
return sum2/s.length();
}
#endif